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1 Introduction 


The Galerkin B-spline method has a number of useful properties (Moser et al. [3]): (i) Thanks 
to the high degree of continuity of the B-splines, the scheme has high resolution similar 
to compact finite-differences (Lele [2]) without the need to separately construct numerical 
boundary schemes, (ii) Constraints such as incompressibility, and coordinate as well as 
flow singularities can be satisfied a priori. Local support allows local constraints as well 
as boundary conditions to be implemented with ease, (iii) Discretized quantities such as 
momentum as well as quadratic invariants such as energy (in the inviscid limit) are conserved, 
(iv) Non-linear terms can be computed without aliasing error, (v) The desired order of 
accuracy is a parameter. 

This work develops another useful property, namely, the ability to do mesh embedding. 
When a gradient becomes large in a certain direction, structured meshes allow one to cluster 
grid lines. However, this is efficient only if the region of large gradient spans the domain. 
This is true only for very simple flows such as laminar shear or shock layers. In more com- 
plicated flows considerable savings can be obtained by using the semi-structured approach 
of embedded meshes of the type depicted in Figure 2. This paper shows how this capability 
is achieved with B-splines. After providing a brief background on one-dimensional B-splines, 
it presents an algorithm for choosing an appropriate set of two-dimensional functions and 
then presents test cases for linear scalar operators. Refinement ratios between mesh blocks 
are arbitrary and the properties listed above are retained, in particular, conservation with 
arbitrary order of accuracy. 1 This order of accuracy is uniform everywhere. In contrast, for 
finite-difference or finite volume methods “simultaneous achievement of both conservation 
and accuracy is very difficult and even impossible in most cases” according to Kallinderis 
(1992). In many such methods the order of accuracy drops at mesh interfaces. Such schemes 
often update each zone separately and interpolate zone boundary information in a separate 
step. This requires down-wind differencing at some interfaces which can be destabilizing. 
Such a procedure is not needed here. 

Since the cost of computing various matrices is incurred every time the mesh changes, 
the present method is not recommended for applications requiring frequent adaptive re- 
meshing. The application we have in mind, namely statistically stationary turbulence, should 
not require frequent re-meshing. The extension to three-dimensional embedding should be 
possible in principle, ignoring the cost of mass-matrix inversion. 

2 Background on one- dimensional B-splines 

For the purposes of this work, a one-dimensional spline is defined to be a polynomial of 
degree d in each interval with d — 1 continuous derivatives across interval boundaries. The 
boundaries of the intervals are called knot-points: the dth. derivative of the spline has a jump 
at the knot points interior to the domain. 

A B-spline is simply defined as a spline which has support over a minimum number 

1 For uniform knots, resolution (property (i)) is tested using operators whose eigenfunctions oscillate 
uniformly everywhere. This is not appropriate for embedded meshes and better tests, perhaps involving 
non-constant coefficients, are needed. 
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Figure 1: Plot of one-dimensional quadratic B-splines on the set of knot points indicated by 
X. All functions, except those near the boundary, have support over three intervals. Note: 
Line types are re-used for different functions. 


of intervals and which is normalized. By equating the number of known continuity and 
normalization conditions to the number of unknown coefficients one finds that the number 
of minimum intervals is d + 1. This is enough information to determine the B-splines for 
all consecutive d + 2 points. Thus a quadratic B-spline has support over three intervals. 
Figure 1 shows the set of quadratic B-splines for the knot points indicated by x . Near the 
boundary, the number of continuity conditions available drops and so does the number of 
intervals of support. To calculate the B-splines ones does not have to actually solve any 
continuity conditions. Rather, one uses a recurrence relation given in de Boor [1] (Chapter 
10) to build up the functions from piecewise constant functions. 

Above, the B-splines were merely defined as splines with minimum support. What makes 
them useful for solving partial differential equations is the result, due to Curry and Schoen- 
berg (see de Boor [1]), that they form a basis for spline functions with the given knots. 


3 Function selection algorithm 

3.1 Mesh Definition 

It is assumed that the computational domain is mapped to a rectangle in £, 77. The case 
of a more general polygon with right angles (such as a backward-facing step) can perhaps 
be treated along very similar lines but this is not presently allowed. For convenience, and 
without loss of generality, the user is required to specify the mesh in terms of sets of points 
which are swept across certain intervals. For instance in Figure 2, the horizontal lines of the 
mesh can be generated by sweeping the three sets of 77 points indicated by • across the £ 
intervals indicated. Similarly, the vertical lines can be generated by vertically sweeping the 
two sets of £ points indicated by x . Just as in the ID case where the interval boundaries are 
knot points, we want mesh lines in the 2D case to represent knot lines, i.e., the lines normal 
to which the selected functions (of degree d, say) have a jump in the dth normal derivative. 
Hence we regard each set of points used to sweep out the mesh as being a knot-set with an 
associated set of one-dimensional functions. The first £ knot set in Figure 2 is the same as 
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Figure 2: Sketch to illustrate the function selection algorithm. All functions are quadratic 
B-splines and span three knot intervals. 


the knot set shown in Figure 1 and has the same associated functions. 

If every £ knot set either contains or is contained by every other £ knot set (and similarly 
for 77), we will say that the mesh is hierarchical. Figure 2 shows a non-hierarchical mesh and 
the same algorithm applies for both types of meshes. 

It is assumed that the mesh thus specified by sweeping points produces closed cells. 
Suppose that a vertically swept point £1 “hangs” at 771, i.e. fails to continue into the next 
sweep. Then to ensure a closed cell 771 should be a knot point in the 77 knot set which sweeps 
over £1 or begins or ends a sweep there. 

As a preliminary, break up the mesh into a set of blocks on each of which the mesh is 
regular. Figure 2 has five blocks. This decomposition is not always unique and should be 
performed to minimize block boundaries. In the present implementation a provisional set of 
blocks is first created on the basis of the sweep intervals. These blocks are then merged with 
neighbors to create larger blocks. 

Below, an algorithm with variants is presented for selecting functions. The first one, 
the intersection procedure, does not allow functions that create knot lines not defined by 
the mesh. For non-integer refinement ratios, this restriction would produce coarse functions 
along block boundaries and the less constrained procedures described in §3.3 should be used. 
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3.2 Intersection procedure 

We want to choose a set of functions B n (£,r)),n = 1,2, ...TV that can represent piecewise 
polynomials of degree d having knot lines that coincide with the mesh. We have no proof 
that the procedure provides a complete basis. Each B n (£,Tj) is constructed as a product 
of one-dimensional B-splines, f(£)g(ri), say. One may consider implementing a brute-force 
algorithm in which all possible pairs of £ and 77 functions defined by the swept knot sets are 
tried and those that create undesired knot lines are rejected. This procedure would not only 
be too costly but does not allow the choice of functions not defined by the given knot sets 
which, we shall see, come into play for non-hierarchical meshes. 

Two-dimensional functions confined to each block are chosen a priori and it remains to 
select functions that penetrate multiple blocks. They are referred to as “spilling functions”. 
The procedure is to consider every function, say /(£), on each sweep of the £ knot sets 
and find suitable functions of 77 as multipliers. The same procedure is repeated for every 77 
function on each sweep of the 77 knot sets. Constant reference is made to Figure 2. 

(1) First, one obtains the 77 knot set from which multipliers for /(£) will be chosen. The 
^-support of /(£) will penetrate a certain range of sweep intervals of the 77 knot sets. 
For instance the function A penetrates sweeps 1 and 2 of the 7/ knot sets. To prevent 
creation of new knot lines, we must choose multipliers from the intersection of the 77 
knot sets associated with the range of sweeps. We will call this the compatible knot 
set of /(£). For instance in Figure 2, the function B (from the first 77 knot set) is 
not a compatible multiplier for A because the resulting product has additional knot 
lines indicated by the dashes. Function C (from the intersection of 77 knot sets 1 and 
2) is a compatible multiplier. Note that if the knot-sets involved in the intersection 
operation are hierarchical (i.e. each set either contains or is contained in another) then 
the intersection operation just chooses the coarsest of the sets. In general, however, the 
compatible function could turn out to be one that does not belong to any of the knot 
sets used in the mesh definition. 

A reviewer suggested the following for clarification. Note that A is the left-most function 
in £ knot-set number 1 that will multiply coarse functions in 77. Thus functions (such as 
A x C) that have the y knot spacing of the coarse block penetrate two £ intervals of the 
fine block. In general the number of intervals of penetration is the polynomial degree, d. 
This is analogous to the fact that in one-dimension, functions with non-zero value and 
derivative at a boundary penetrate d intervals into the interior of the domain. 

(2) Only those functions in the compatible knot set that have support in the sweep interval 
of the knot set of /(£) are relevant; we refer to these as compatible functions. 

(3) Next, further limit the compatible functions in order to prevent block-confined functions, 
which have been chosen a priori , from being selected. The support of /(£) x its 77 sweep 
interval is either confined to a block or it is not. 

(a) If it is confined to a block then only those compatible functions that cross the block 
boundaries in 77 are allowed. 

(b) If it is not confined to a block (as is the case for function A), then all compatible 
functions are allowed. 

(4) In order to prevent the selection of 2D functions with additional knot lines, it is necessary 
that the compatibility be mutual. For instance D is a compatible function of A, but A 
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is not a compatible functions of D. This is because function D penetrates sweeps 1 and 
2 of the £ knot sets but the knots of function A do not belong to their intersection. In 
particular, the following test is applied: do the knots of /(£) belong to the compatible 
knot-set of g(rj)l 

The mutual compatibility test needs to be relaxed at non-hierarchical corners. Consider, 
for instance, the region near 0 where two fine regions meet at a corner. The only 
functions that have support at the corner 0 that do not result in additional knot lines 
are functions such as E x D, neither of which belongs to any knot-set used in the 
definition of the mesh. Since only functions on the knot sets used in the mesh definition 
seek suitable multipliers, such a product would never be chosen. This is overcome by 
either of two modifications, denoted as Ml and M2. If the mutual compatibility test fails 
for a prospective 2D function containing a block corner and the g(r j) does not belong 
to any of the g knot sets penetrated by the support of /(£) then: (Ml) g(r)) is chosen, 
resulting in the creation of new knot lines, or, (M2) The product of g(g) and all the 
functions on the compatible knot set of g(rj) that have support at the corner are chosen. 
In this case no new knot lines are created in the context of the intersection procedure. 
The unmodified procedure is denoted as MO. 

(5) Finally, a check is made that a 2D function selected is unique. 


3.3 Less constrained procedures 

It may be desirable to change resolution gradually, in a non-integer fashion as shown in 
Figure 2 between the left and right halves of the mesh. If compatible functions for F were 
chosen according to the intersection procedure, the very coarse function G would result (its 
knots are indicated by o). To avoid this, one can dispense with the requirement that no 
new knot lines be created. One simple way of doing this is to use the “densest” instead of 
the intersection operation to determine the compatible knot set (i.e. choose the set with the 
most knots). Another is to use the “densest by sweep (DS)” operation, i.e., assemble the 
compatible knot set by taking the densest set within each sweep. In this case a function 
such as H would be a valid multiplier for F and the product would create the additional 
knot lines indicated by dots. The function H is also mutually compatible with F because 
the knots of F belong to the (single) £ knot set penetrated by H . The mutual compatibility 
test was needed in the intersection algorithm to prevent additional knot lines. Here it is 
necessary to prevent selection of multiple types of products in the same region. For instance, 
function C would seek multipliers at some point. One sees that F would not be a mutually 
compatible function because the knots of C do not belong to the DS of the vertical knot sets 
on which F has support. This is rightly so because F X H is the type of product we have 
selected for this region. 

In the tests that follow the algorithm used will be denoted by a prefix: I (intersection), 
D (densest), or DS (densest on a sweep by sweep basis). This will be followed by a suffix 
(MO, Ml or M2) to denote the modification to the mutual compatibility test. 

Algorithm summary: Perhaps the following summary will aid the reader in holding the 
algorithm firmly in mind. For every £ function on the knot sets whose £ support x g sweep 
interval is confined to a block, pick as multipliers only those mutually compatible g functions 
which cross the block boundaries in g since the rest produce 2D functions confined to the 
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block and have been chosen a priori. For a ( function whose £ support x r/ sweep interval 
is not confined to a block, choose as multipliers all mutually compatible functions. Repeat 
the procedure for all rj functions defined by the given knot sets. If required, relax the test of 
mutual compatibility according to Ml or M2. 

3.4 Matrix structure of linear operators 

In a Galerkin (weighted residual) method, linear operators (such as the Laplacian and ad- 
vection operator in the present examples) give rise to matrices of the following general form: 

(CiB m ,C 2 B n ) = J C 1 B m C 2 B n dxdy, m,n=l,2...N, (1) 

where C\ and £ 2 represent linear differential operators. Similarly, operators with a quadratic 
non-linearity give rise to integrals of triple products. All matrices need be computed once. 
Since differentiation does not alter the support region of a function, all linear operators 
produce matrices with the same structure that depends on which pairs of functions overlap. 
The template for this structure is constructed once. Pairs of functions confined to each 
block overlap in a simple way and produce the structured part of the matrix: it has (d + l) 2 
diagonals for a symmetric operator and (2d + l) 2 diagonals for a non-symmetric operator. 
Spilling and block-confined functions that overlap produce scattered matrix elements as do 
two overlapping spilling functions. With the template in hand, matrix elements are computed 
using Gauss quadrature with enough points to ensure exact integrals. 

Matrix equations were solved using the preconditioned congugate gradient method which 
requires that the left hand side matrix multiply a vector in each iteration. The action of 
linear operators is also a matrix multiplication hence it is important to make this efficient. 
The multiplication of each diagonal vectorizes over the row direction while the template 
for storing the scattered elements ensures that their multiplication can be accomplished in 
a certain number of vector multiplies. This number is the maximum number of scattered 
elements in a row. The template is rearranged to avoid the most obvious bank conflict, 
namely that of accessing the same element of the multiplying or resultant vector within a 
certain number of CPU clocks. 

4 Test cases 

First consider scalar advection at 45° to the x-axis: 

u t t + cu tX > = u )t + c(u yX + u tV ) = 0 on the unit square, (2) 

where x' is a coordinate at 45° to the x-axis. The initial profile u(x', 0) is a Gaussian pulse of 
unit height and e -1 half- width of .15. The advection speed, c, is set to unity. The boundary 
condition imposed on u it at the in-flow (left and bottom) boundaries corresponds to uniform 
propagation of a Gaussian pulse. Henceforth, we use as weight functions those B-splines, say 
B m (x, y),m = 1,2,... N 0 , which vanish where boundary conditions are applied and represent 
the solution as 

No N 

u(x,y,t) = Y^a n (t)B n (x,y)+ ^ a n (t)B n (x,y). (3) 

n=l n=iVo+l 
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No. of coarse intervals 


Figure 3: (a) Local error for the advection equation after the wave has propagated a distance 
of 0.6. Contour min., max. and inc. = (—.0065, .0060, .0005). The error is sampled on a 150 
x 150 grid, (b) Convergence test for the advection equation. Same instant, mesh type and 

sample points as (a). , Maximum error (Loo norm); , Average of absolute error 

(Li norm); , reference lines with slope of -3 and —4; A , for quadratic splines; ■ , for 

cubic splines. 


Denoting the inner product as (., .), the weak formulation reads: 

N 0 N N 

} B n )a n = c ^ ](B m , B nx 4* B n v )a n S . m = 1, 2, . . . iV 0 (4) 

n=l n=l n=JV 0 +l 

where the coefficients in the last term are known from projecting the boundary specification 
of u tt . In order to make time integration errors smaller than spatial discretization errors, 
this equation is advanced with a sixth-order Runge-Kutta scheme with cfl = cAt/Ax^n = 
0.5. Figure 3a shows the local error obtained using cubic splines and algorithm DM0 for 
a non-integer refinement ratio of 5/4. At the instant shown, the pulse has propagated 4 
half-widths and the peak of the pulse lies on the diagonal of the fine block. The error 
is perfectly symmetric about the 45° line, smooth, without any peculiarities at the mesh 
interface, and attains a maximum of 0.65% even on the rather coarse mesh. Figure 36 is the 
result of a convergence test using quadratic and cubic splines. Aside from a little curvature 
and flattening for the cubic case, the approximation-theoretic convergence rate (polynomial 
degree +1) is obtained. 

Next consider the Poisson equation for the streamfunction given the vorticity, u>: 

V 2 if)(x,y) = —u>(x,y) on the unit square ft with if) = g(s) on dfl. (5) 

In order to eliminate unknown boundary terms in the weak formulation and to impose the 
Dirichlet boundary condition strongly we again use as weight functions all the B-splines 
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No. of coarse intervals 


Figure 4: (a) Solution to the three-vortex Poisson equation sampled on a 100 x 100 grid. The 

actual mesh is shown. Computed solution: , negative contours; , positive contours. 

Exact solution: , negative contours; , positive contours., (b) Convergence for the 

three-vortex Poisson equation test. Error sampled at 512 x 512 points: □ , Maximum 

error for the embedded mesh; + , Maximum error for a uniform mesh having the finest 

spacing of the embedded mesh. Error sampled at all the “mesh-points”: A , maximum 

error (L^ norm); o , r.m.s. error ( L 2 norm); • , Average of absolute error (Li 

norm). , reference lines with slopes of —3 and —4. 


which vanish at the boundary. The weak formulation reads: 

No N 

VB n ) = V£ n ), m = l,2,...AT 0 (6) 

n=l n=No + 1 

The test vorticity field in the present example consists of three axisymmetric Gaussians of 
alternating sign and graded intensity superposed with images to make the left and bottom 
boundaries impermeable walls. The exact streamfunction due to each vortex can be expressed 
in terms of the exponential integral. In the numerical solution the condition tj) = 0 is 
imposed on the left and bottom walls and the exact solution is imposed on the top and right 
boundaries. 

Figure 4 a shows that even for the very coarse mesh, the agreement between the exact 
and numerical solutions is excellent (the intersection algorithm (IMO) with quadratic B- 
splines was used). For comparison, the solution was also obtained for a uniform mesh having 
everywhere the mesh spacing of the finest block of the embedded mesh. The maximum error 
occurs in the intense vortex near the corner and it has virtually the same value for the two 
meshes (compare o and + in Figure 46; the error was sampled on a 512 x 512 set of points 
for both meshes). 

The rest of the curves in Figure 46 show that when the error is evaluated at the “mesh- 
points” (i.e. the intersection points of the knot lines), two of the norms converge at fourth 


8 





order which is one order higher than the approximation theoretic result. This phenomenon 
is called super-convergence and deserves a brief comment. Thomee [4] showed that in one- 
dimension and with periodic boundary conditions, the Galerkin B-spline method exhibits 
super-convergence at the knot points. In particular for the heat equation the convergence rate 
is 2d while for the advection equation it is 2d + 2. Exact time integration is assumed in both 
cases. Various numerical tests were performed and they indicated that super-convergence 
was weak for two dimensions with non-periodic boundary conditions: Figure 4 6 represents, 
indeed, an exceptional case. For the advection equation there is no super-convergence even 
on a uniform mesh. For the Poisson equation on a uniform mesh, the odd-degree splines 
considered did not exhibit superconvergence, while among the even-degree splines only the 
quadratic functions displayed a rate consistent with Thomee’s result. Quartics converged 
at sixth order compared with fifth order for approximation theory and the eighth order of 
Thomee’s result. 

In the advection test the local error had no peculiarities near the mesh interface. The 
same is not true for the Poisson test. Figure 5a shows that the maximum error occurs in 
the region of the intense vortex nearest the corner, as one would expect, but increases again 
in the form of positive and negative layers along the interface with diminishing amplitude 
normal to the interface. A regular mesh (without embedding; see Figure 56) but with the 
same change in spacing normal to the interface as the embedded mesh gives similar features 
in the error but the peak is a little smaller (70% of the value in the embedded mesh). In 
Figure 5c normal spacing is kept uniform while tangential resolution changes. The pattern 
of the error near the interface resembles swords pointing normal to the interface with the 
positive peak being half the value in the original embedded mesh. Figure 5 d shows that 
near the interface of interest, the two errors very nearly add to give the error in the original 
embedded mesh. 

The final test is one of robustness: we consider a problem for which a uniform mesh 
is optimal and study how much the resolution degrades when the mesh is dislocated (as 
shown in Figure 66) and the embedding procedure is applied. A uniform mesh is optimal 
for problems in which the solution oscillates uniformly everywhere in the domain. The 
eigenvalue problem for the Laplacian operator is one such problem: 


O / 

VV = Xxp on a rectangle of unit width and height h , with — — = 0 on the boundary. (7) 

on 


The von Neumann boundary condition was chosen because it makes the boundary term of 
the weak formulation vanish. The square, h 2 , of the height of the domain is chosen to be 
irrational ( h 2 = y/2) to avoid degenerate eigenvalues. The exact eigensolutions are: 


A m ,n = -7r 2 (m 2 + (rc/ft) 2 ), = cos(m7ra:) cos(n7ry), (8) 

The extent to which a numerical method is able reproduce the exact eigenvalues and eigen- 
functions is a test of its resolution of the Laplace operator (with the given boundary condi- 
tions) across all scales. We compare performance on a uniform 10 x 10 interval mesh (without 
embedding) with performance on the dislocated mesh shown in Figure 66. 

The pairing of a numerical eigensolution with an exact one is somewhat arbitrary for the 
inaccurate eigensolutions. The procedure used was to evaluate each numerical eigenfunction 
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Figure 5: Local error on three different meshes for the three-vortex Poisson equation test. 
In all plots, contour min., max., inc. = (—.00475, .00425, .00050). This makes the lowest 

contour levels ±.00025 rather than zero. , positive values; , negative values. The 

error is evaluated on a 300 x 300 grid, (a) Mesh embedding with abrupt changes in spacing 
normal and tangential to the mesh interfaces, (b) A non-embedded but non-uniform mesh 
with abrupt changes in normal resolution only, (c) Mesh embedding with an abrupt change 
of spacing in the tangential direction in the region of interest, (d) The result of adding (b) 
and (c). 


Table 1: Number of eigenvalues (n\) and eigenfunctions (n„) with error less than the specified 
tolerances. Relative error is used for eigenvalues; Li error is used for eigenfunctions. 


Mesh type 

10% error 

tolerance 

1% error 

tolerance 

No. of degrees of freedom 

Uniform 

n\ = 90 

n v = 63 

o 

CO 

II 

£ 

n v = 22 

144 

Dislocated 

n\ = 93 

n v = 57 

n\ = 29 

CM 

CM 

II 

> 

c 

151 
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Figure 6: Local error for the eigenfunction for m = n = 6. In both plots contour min., max., 
inc. = (—.19, .17, .04) which makes the lowest contour levels —.03 and .01 rather than zero. 

, positive values; , negative values. The error is evaluated on a 50 x 50 grid, (a) 

Uniform mesh (b) Dislocated mesh. 


on a 50 x 50 grid, normalize it to have a value of unity at the origin and pair it with the 
closest (in the discrete L-2 sense) unpaired exact eigensolution. The pairing proceeded in 
order of increasing numerical eigenvalue. To assess resolution, the number of eigenvalues and 
eigenvectors which satisfy a given error tolerance for the two meshes are shown in Table 1. 
For a tolerance of .01, the number of “good” eigensolutions is only slightly lower for the 
dislocated mesh. The overhead of extra functions at the dislocation increases the number 
of degrees of freedom from 144 to 151. Therefore the ratio of good eigenvectors to the total 
number of degrees of freedom degrades by 5% for a tolerance of .01. For a tolerance of .10 
this ratio degrades by 14%. Figure 6 compares the error in the eigenfunction which has 
three wavelengths over the ten intervals in each direction. The dislocated mesh actually has 
a maximum error which is smaller by 3% but L 2 error which is larger by 3%. 


5 Conclusion 

A technique has been developed for achieving two-dimensional mesh embedding in the con- 
text of a Galerkin method with B-splines as basis functions. The results of test cases are 
encouraging and work is under way to apply the technique to the incompressible Navier- 
Stokes equations for three-dimensional flow in two-dimensional curvilinear coordinates. 
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